Paired chiral spin liquid with a Fermi surface in S — 1 model on the triangular lattice 
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Motivated by recent experiments on Ba3NiSb2 0g, we investigate possible quantum spin liquid 
ground states for spin 5 = 1 Heisenberg models on the triangular lattice. We use Variational 
Monte Carlo techniques to calculate the energies of microscopic spin liquid wave functions where 
spin is represented by three flavors of fermionic spinon operators. These energies are compared 
with the energies of various competing three-sublattice ordered states. Our approach shows that 
the antiferromagnetic Heisenberg model with biquadratic term and single-ion anisotropy does not 
have a low-temperature spin liquid phase. However, for an S(7(3)-invariant model with sufficiently 
strong ring-exchange terms, we find a paired chiral quantum spin liquid with a Fermi surface of 
deconfined spinons that is stable against all types of ordering patterns we considered. We discuss 
the physics of this exotic spin liquid state in relation with the recent experiment and suggest new 
ways to test this scenario. 



PACS numbers: 71.27. 



75.10.Jm, 75.10.Kt, 75.30.Kz 



-f— > 



i 

C 

O 

o 



m 

(N 
rn 

O 



I. INTRODUCTION 

Quantum spin liquids (QSL) are interesting states 
of matter with long-range entanglement that may ex- 
hibit exotic properties like unbroken lattices symme- 
tries at low temperature, quasiparticle fractionalization, 
emergent gauge fields, braid statistics, and chiral edge 
modes. 1 '- The existence of QSL or resonating- valence- 
bond (RVB) states in two dimensions was first conjec- 
tured by Anderson as possible low-temperature phases 
of the spin-1/2 antiferromagnetic Heisenberg model on 
the triangular latticed Shortly after, a fascinating re- 
lation of RVB states with high-temperature supercon- 
ductivity was uncovered: Upon doping, some quan- 
tum spin liquids are expected to give rise to uncon- 
ventional superconductivityiir£ So far, several experi- 
ments found indication of spin liquid behavior in a num- 
ber of geometrically frustrated two-dimensional spin-1/2 
antiferromagnets.— ~— Spin systems with higher values of 
spin, however, usually show a strong tendency towards 
long-range ordering and lattice-symmetry breaking at 
low temperature. 

Last year, highly surprising experimental results 11 
found spin-liquid behavior in new structural phases of 
Ba3NiSb20g. In the so-called 6H-B phase, obtained 
through a high-pressure treatment of this antiferromag- 
netic insulator, the Ni +2 -ions, carrying effective spin 
5 = 1, arrange in presumably weakly coupled layers of 
triangular lattices. No magnetic ordering was observed 
down to 0.35 K despite a large Curie- Weiss temperature 
of &cw — —75 K, and the magnetic susceptibility (af- 
ter substraction of orphan spin contribution) was found 
to saturate at low temperature T. Furthermore, mea- 
surement of the magnetic specific heat found Cm oc T. 
These properties are highly unusual for an insulator but 
are typical for metallic states. For example, spin-wave 
theory for conventional long-range ordered states predict 
a specific heat Cm oc T 2 in two space dimensions. 12 



So far, a number of theoretical attempts have been 
made to explain these experiments. Two possible spin- 
liquid candidates were proposed In Ref. 14], a de- 
composition of spin 5 = 1 operators in terms of four 
flavors of fermionic spinons and their possible mean-field 
states were conjectured. Such a decomposition is most 
natural in the case of a two-orbital Hubbard model with 
not too strong interactions (Hund coupling) between the 
electrons. The minimal number of spinons required to 
represent spin 5 = 1 is t/treei 15 ' 16 On the basis of this 
three-fermion representation, an exotic QSL state was 
proposed by some of us in Ref. [l3| that well reproduces 
the phenomenology of the experiment on Ba3NiSb2 0g. 
A less exotic scenario was recently proposed in Ref. [17| , 
where strong inter-layer couplings between the Ni 2+ spins 
drive the system to a quantum critical point. These au- 
thors predicted a T-linear specific heat in some temper- 
ature range. 

Currently, the details of the effective spin model de- 
scribing Ba3NiSb2 0g are not known. In this paper, we 
will not propose a realistic effective spin model for this 
material. Instead, we want to investigate two families of 
promising antiferromagnetic triangular-lattice spin-one 
models at the variational level. The aim is to deter- 
mine whether, variationally, the natural quantum spin- 
liquid candidates have a chance to win over long-range 
ordered ground states in these models. First, we consider 
the bilinear-biquadratic Heisenberg model with single-ion 
anisotropy term. In this model, we do not find evidence 
for a low-temperature QSL phase. We further propose 
an SU(3) symmetric model with three-site ring-exchange 
terms. In this model, for strong ring-exchange terms, we 
find that an exotic spin liquid state is stabilized. We dis- 
cuss the phenomenology of this state and propose further 
experimental tests of this scenario. While the theory for 
5 = 1/2 QSL is well developed and has a long history, 
much less is known about spin liquids for 5 = 1. Here 
we present new methods and results on this problem^ 
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This paper is organized as follows. In the next section, 
we introduce the representation of spin S — 1 in terms 
of three fermionic spinon operators. Section ITTll describes 
all spin liquid wave functions as well as the long-range 
ordered states that we are considering. In Sec. IIV1 we 
introduce two spin models and present the variational 
results we found for these models. In Sec. [V] we discuss 
the chiral edge modes corresponding to the chiral QSL 
state that we found to be stabilized in one of the mod- 
els. Section IVT1 discusses the response function and other 
physical properties of this state and, finally, we conclude 
in Sec. En! 



II. SPINON REPRESENTATION 

To construct spin liquid states for spin S — 1, we fol- 
low an approach similar to the one outlined in Ref. [l5j . 
We write the spin operators in terms of three flavors of 
fermionic spinons, f a , in the following way: 



&a — i^abcflfc j 



(1) 



where a € {x, y, z}. In this paper, repeated indices are 
always summed over. We choose to work with operators 
f a that create spin states |o) in the time-reversal invari- 
ant basis, i.e., 



\x) = ^=(\i) + \i)), 
|tf> = ^=(|i>-|I», 

|*)=-*|0), 



(2) 



where |1), |1), and |0) are 5 2 -eigenstates with eigenvalues 
±1 and 0, respectively. 

By representing spin in terms of fermions, we have en- 
larged the Hilbert space. The fermion operators act in 
the eight-dimensional Fock space while the original spin 
space is three dimensional. In order to recover the physi- 
cal subspace, a local constraint on the fermionic occupa- 
tion number has to be enforced, 



In the time- reversal invariant basis ([2]) , the quadrupo- 
lar operators, defined as Q a b = (S a Sb + S S a — 3)/2, 
acquire a particularly simple form.— In the particle rep- 
resentation (Nf — 1), we have 



S a Sb — S a b ~ flfb ■ 



(6) 



• — y ] fa fa 



N 



f 



(3) 



In order to analyze a particular spin S = 1 lattice 
model within the spinon representation (JTJ) , one may 
start by decoupling the spinon-interaction terms with the 
help of a Hubbard-Stratonovich transformation. To im- 
plement the constraint ([3]) and the symmetry properties 
(HJ and iJSJ in this theory, a compact gauge potential for 
the local symmetry group has to be introduced in the 
path integrali&Si This procedure enables derivation of a 
low-energy effective theory for possible spin liquid phases 
but does not address microscopic stability for a particu- 
lar Hamiltonian. Better suited for this purpose is a vari- 
ational wave function approach. In this approach, un- 
physical states are removed by hand from wave functions 
that correspond to possible low-temperature phases of 
the theory. This allows the construction of a new class of 
genuine microscopic variational wave functions for spin- 
one models. Determining the best variational state for a 
spin model then provides guiding information about the 
low-temperature phase of the model. In the present pa- 
per, we first follow this approach, and discuss possible 
microscopic Hamiltonians. We also use the effective field 
theory to discuss some properties of the proposed QSL 
states. 



III. VARIATIONAL WAVE FUNCTIONS 

In this section, we introduce two classes of microscopic 
variational wave functions for spin S — 1 on the triangu- 
lar lattice. First, we describe quantum spin liquid wave 
functions that do not break the space group symmetries 
of the lattice. Second, we outline a general approach for 
constructing competitive long-range ordered states that 
have an enlarged unit cell. 



A. Quantum spin liquid wave functions 



Both particle (Nf = 1) or hole (Nf = 2) subspaces can be 
chosen. Furthermore, the spin operator remains invariant 
under the transformations 



fa 



i<t> 



fa 



and 



fa f a 



(4) 



(5) 



Eq. ([5]) is a particle-hole transformation and the con- 
straint ([3]) is changed according to Nf M> 3 — Nf. Hence, 
the local symmetry group for this representation of spin 
5 = 1 operators is the semi-direct product U(l) x Z2. 15 



We start by writing down quadratic "trial" Hamiltoni- 
ans in terms of the spinon operators, 



H qsl = ^{ S / a y Qj + + h.C.} 

(i,3) 

~ Ma fljfaj ■ 
3 



(7) 



The sum goes over the nearest-neighbor links of the 
triangular lattice. In this trial Hamiltonian, the emergent 
gauge fields that would be present in the corresponding 
low-energy theory are omitted. Particular values for the 
mean-field parameters s = ±1, A" fc , and jjL a represent 
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possible low-temperature QSL phases, s = ±1 corre- 
sponds to flux of 7r or zero through all triangles of the 
lattice. Next, we assume unbroken spin-rotation sym- 
metry around the z-axis, and we focus on iS* ot eigen- 
states with 5^ ot = 0. Note that under spin rotations, 
fj = {fx, fy, fz)j transform as real vectors [i.e., fj trans- 
form in the adjoint representation of SU{2)]. Further- 
more, we restrict ourselves to states that have a single 
site per unit cell and that do not break the space group 
symmetries of the lattice (translations, rotations, and in- 
version). In this situation, the following cases exhaust 
the possible QSL candidates on the triangular lattice: 

(i) [7(1) state: = 0. 

(ii) Equal-flavor pairing: Af? ^ 0, Af/ = Af/ + 1 0, 
A" h = otherwise. 

(iii) x-y pairing: A?/ = — Af/ ^ 0, Afj = otherwise. 

The chemical potentials for x- and y-fermions are cho- 
sen to be identical, [i x = fi y . Other possible pairings 
Afj than the ones considered in (ii) or (iii) violate our 
symmetry requirements^ 

On the one hand, for Af/ = Af/ = Af?, equal- flavor 
pairing (ii) corresponds to spin-one singlet pairing. The 
pairing term in (JT]) creates a state fj ■ fj\0) that is in- 
variant under all spin rotations; hence, it is a singlet. In 
general, for Af* — Af/ ^ Afj, the state is not an eigen- 
state of (S|f ) 2 = {Si + Sj) 2 . However, for Af/ = Af/ = 
— Af//2 one can check that (S 1 */') 2 = 6; therefore, this 
bond operator creates a spin-one quintuplet. On the other 
hand, the x-y pairing bond operator Ay/^/y, (iii), cre- 
ates a spin-one triplet. To see this, let us denote the 
state by = {\xy) i} - \yx) lj ) oc (|ll)y - |ll)y)- Since 
(S\f) 2 =4 + 2Si- Sj, and [Si ■ Sj + 1]|% = 0, we have 
(S^) 2 = 2. Note that the total spin per site for all these 
QSL states is small in the thermodynamic limit. We have 
^((5' tot ) 2 )/A r ~ 1/VN where N is the number of sites 
and S toi V 

Due to the anti-commuting spinon operators, the pair- 
ing parameters Afj must have particular symmetry prop- 
erties under inversion of the link direction For 
equal-flavor pairing (ii), we have Afj = —A"", i.e., the 
pairing is odd under space inversion. For x-y pairing (iii) , 
we have Af/ = A x f , i.e., the pairing is even under space 
inversion. This is in contrast to S = 1/2 spin liquids, 
where singlet pairing is even while triplet pairing is odd 
under space inversion. 

In order to obtain a microscopic variational QSL wave 
function, we take the ground state |^o) of ([7]) and apply 
the Gutzwiller projector Pcin-j = 1), enforcing single 
occupancy on each site and thereby removing unphys- 
ical components. In this way we construct a genuine 
spin-one resonating- valence-bond (RVB) spin- liquid wave 
function, generalizing similar approaches to S — 1/2 spin 
liquids^ In this paper, we choose to work in the micro- 
canonical formalism where the fermion number is held 



fixed, i.e., we project the wave function to a fixed total 
number of spinon flavors, N a — Y]j n a j, 

\N)=P N P G (nj- = l)|ifo), (8) 

with N = (N x ,N y ,N z ). Since N x = N y (to maintain 
spin-rotation symmetry around the z-axis) and from the 
local constraint we have 2^ + N z = N , where TV is the 
number of lattice sites (N = 12 x 12 in most of our calcu- 
lations). Expectation values in RVB wave functions (|8]) 
can be calculated numerically within Variational Monte 
Carlo (VMC) techniques^ More technical details on our 
numerical scheme are given in the appendices. 

The possible complex phases (pairing symmetries) of 
Afj in (ii) and (iii) are restricted by the rotation sym- 
metries of the lattice: Let us denote the nearest-neighbor 
links of the triangular lattice by 1 = (1,0), and 2,3 = 
(±1, v / 3)/2. For equal- flavor pairing (ii), the pairing sym- 
metry can be real /-wave with A? a = — A| a = A| a , or 
complex p x + ipy-wave (p+ip), with A? a = A? a e~ l7T / 3 = 
A| a e~ i27r / 3 . For x-y pairing (iii), the possible pairing 
symmetries are extended s-wave with A?" = A| y = A^ y 
and d x + id y -wsve (d+id) with A^ = A^ y e "' i27r / 3 = 

A x y e~ l47r / 3 . Higher angular momenta would require 
spinon pairing between further-neighbor sites, which we 
choose to exclude from the present study^ 

Symmetry of the QSL states © under lattice rotations 
forbids mixing of different types of pairing symmetries in 
the Hamiltonian ([7]). For example, a state where the f z 
spinon is paired with /-wave, and f Xl f y are paired with 
p+ip pairing symmetry breaks lattice rotation. Similarly, 
in the x-y paired QSL (iii) , f z must remain unpaired un- 
less lattice rotation symmetries are broken. The reason 
is the following: After performing a lattice rotation on 
the mean-field Hamiltonian ([7]), one would like to find a 
gauge transformation ((4]) that brings it back to the origi- 
nal form. If such a gauge transformation exists, then the 
corresponding spin wave function ([8]) is unchanged by the 
rotation (after Gutzwiller projection and up to a phase). 
However, since all three spinons transform with the same 
U(l) phase, such a gauge transformation can only exist 
when all spinons have identical pairing symmetries! 24 ' 25 

The QSL states have the following properties: Ex- 
tended s-wave and /-wave states respect parity P (reflec- 
tion on a symmetry axis of the lattice) and time reversal 
symmetry O. The p+ip and the d+id states, however, 
break both P and O, but conserve the product OP. In 
this sense, they can be termed chiral spin liquids^ al- 
beit for spin S = 1. The p+ip state is fully gapped and, 
therefore, a conventional topological state of matter. The 
d+id state, on the other hand, represents a new class of 
paired chiral states in two dimensions that exhibit both 
P- and ©-symmetry breaking and a gapless bulk Fermi 
surface at the same time. These exotic properties will be 
discussed in more detail below and in later sections. 

In the U(l) spin liquid (A afc = 0), all three spinon 
excitations have a Fermi surface. This corresponds to 
the Coulomb phase of the emergent U{\) gauge theory 
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where the photons are massless. The paired states with 
A afc 7^ correspond to "Higgs" phases where the global 
U(l) symmetry is spontaneously broken and the pho- 
ton acquires a mass^i Among the equal-flavor pairing 
states, the /-wave state has gapless nodal points in the 
spectrum while the p+ip QSL is fully gapped. In the x-y 
paired QSL, the spin excitations are gapped. However, 
the nematic (S z = 0) excitations form a gapless Fermi 
surface of weakly interacting (and therefore deconfined) 
spinons. Specific heat and spin-susceptibility of an x-y 
paired (triplet) QSL are consistent with the recent ex- 
periments on Ba3NiSb2 0g. 13 

The variational parameters we are using for the micro- 
scopic QSL wave functions are the amplitudes |A ah | for 
all pairing symmetries discussed above, and the chemi- 
cal potentials \i x and \i z . Furthermore, we consider the 
cases s = ±1 in 0, corresponding to the presence or the 
absence of 7r-flux through the triangles of the lattice. For 
the paired states, N z is used as an additional variational 
parameter (independent of \x z ; see Appendix [Bl for more 
details). 



B. Long-range ordered states 

In order to make reliable statements about the low- 
temperature phase of a spin model, the energies of QSL 
wave functions have to be compared with competitive 
long-range ordered states. Here, we consider natural 
ordering patterns that are suggested within a simple 
product-state ansatz (e.g., a 120° magnetic ordering in 
the case of the antiferromagnetic Heisenberg model on 
the triangular lattice). The QSL wave functions (JSJ are 
highly correlated states. To be able to compare the varia- 
tional energies, we also need to introduce nontrivial quan- 
tum correlations to the ordered states. 

Here, we use the following two complementary schemes 
to introduce quantum corrections on top of long-range or- 
dered product states. The first approach builds on the 
fcrmionic representation and gauge theory description of 
the spin model. Long-range ordered phases can be cap- 
tured within the following quadratic trial Hamiltonian, 

fljfaj ■ 

3 

Similar to the QSL wave functions ([5]), the Gutzwiller- 
projected ground state of (|9|) serves as a variational state. 
The normalized complex vectors dj specify a particular 
spin-one ordering pattern. The variational parameter h 
interpolates from the U(l) spin liquid (h = 0) to the 
product state \ip p ) = Y\jJ2a d< j\ a )j wnen h — > °°- As 
before, we set fi x = fi y ; fx x — (x z is taken as a variational 
parameter and we consider 7T- and 0-flux states by s = 
±1. 



Another route to constructing correlated long-range 
ordered wave functions is to apply spin Jastrow factors 
to a product state. The analysis of such wave functions 
for the spin- 1/2 antiferromagnetic Heisenberg model on 
the triangular lattice was pioneered by Huse and Elser 
in Ref. [28|. For that model, Huse-Elser wave functions 
were found to give low variational energies, comparable 
to exact energies on small clusters. A generalization of 
the Huse-Elser wave function to the case of spin 5=1 
can be written as 

\J) = exp(- Y,{P{S zi S zj ) + 1 {S zl S ZJ ) 2 }M v ) . (10) 

(i,3) 

Here, \ip P ) is a product state of spin one. In this paper, 
we restrict ourselves to nearest-neighbor Jastrow factors, 
and take /3, 7 to be real variational parameters. 
A general spin-one product state can be written as 

3 a 

where |o) € {\x), \y), \z)} span the local Hilbert space; see 
Eq. Let us write d = u + iv, where u and v are real 
vectors, and consider the one-site state \tp) = J2ad a \ a )- 
We can always take d = (d x , d v ,d z ) to be normalized and 
u • v = (choice of phase) . The spin expectation value 
in this state is given by 

(S) = 2uAv. (12) 

If d is real, the corresponding state is a spin-nematic with 
(S) = 0. In this case, d is called the director and we have 

(S 2 a ) = l-(d a f. (13) 

On the other hand, w = v = 1/2 corresponds to a spin 
coherent state where \ (S) \ = 1 is maximal. 19 

The fermionic states © and the Huse-Elser wave func- 
tions pop are two quite general and complementary ways 
to introduce nontrivial quantum fluctuations on top of 
spin-one product states (fTT|) . Although additional vari- 
ational parameters can be built into the product state 
itself, we need to choose a suitable (family of) product 
states (specified by dj) to start with. As the example by 
Huse and Elser— suggested, good ground-states energies 
can be obtained by choosing the product states that min- 
imize the energy of the spin model. Below we will use a 
similar choice. 

IV. MODELS AND VARIATIONAL RESULTS 

A. Bilinear-biquadratic Heisenberg model with 
single-ion anisotropy 

We start by considering the simplest extension of the 
spin-one Heisenberg antiferromagnet on the triangular 
lattice, 

H KD = Y J {S i ■ Si + K{S t ■ Sj) 2 } + S% , (14) 

(iJ) 3 
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where we set the Heisenberg exchange energy J = 1. 
In this study, we want to restrict ourselves to the pa- 
rameter range \K\ < 1.5 and |D| < 1.5. For D = 0, 
it is known that the ground state of this model ex- 
hibits 120° antiferromagnetic order when K < 1, and 
90° antiferro-nematic (also called antiferro-quadrupolar) 
order when K > l , 29 i 30 For large easy- axis anisotropy 
the ground state is a ferro-nematic product state 
with S Z j = on every site. 

For intermediate values of D, one may expect an x-y 
paired QSL to be stabilized in this model. Since S^j = 
1 — n z j, the single-ion anisotropy D acts as a chemical 
potential for the f z spinon. For non-zero D, the Fermi 
surfaces of f z and of f x , f y in the U(l) state are expected 
to mismatch, and it is conceivable that f x and f y pair 
while leaving f z with a spinon Fermi surface. Indeed, 
unconstrained mean-field theory in Ref. [l~3| found that, 
for K < 0.5, the p+ip QSL wins, while for K > 0.5, the 
d+id state with a spinon Fermi surface is the most stable 
QSL candidate (However, in contrast to above intuition, 
the phase boundary showed only weak dependence on D). 
In the present paper we find that the variational energy 
of ordered states is always lower than the one of the QSL 
states when the local constraint ([3]) is taken into account 
exactly. 

A variational study of the model (fl4l) at the level 
of product states performed in Ref. [3l| suggested a 
three-sublattice orde ring pattern, generalizing the order- 
ing pattern of Ref. [29] to D ^ 0. Motivated by this 
proposal, we choose the following two classes of three- 
sublattice product states as an input to our correlated 
ordered states discussed above. First, we consider an an- 
tiferromagnetic (AFM) state where the spins (V'plSjIV'p) 
have a constant length and lie in a common plane at 
an angle of 120° to each other on nearest-neighbor sites. 
The average spin length, |(?/'p|S'j|'0p)|, is taken as a vari- 
ational parameter. In the notation of Eq. (|llj) . the spin 
states on sublattices A, B, and C of the triangular lattice 
are written as: 

dj e A = (0, —i sin 77, cos 77) , 

V3i i ( 15 ) 

djeB.c = (±— sin 77, -sin 7/, cos 77), 

where 77 € [0,7r/2] is a variational parameter. Using 
Eq. (fT2jl . one may check that this state corresponds to 
120° antiferromagnetic ordering in the x-y plane with 
K^iplS^l^p)! = sin2?y. We also consider the same order- 
ing in the x-z planed For 77 = 7r/4, each site is in a spin- 
coherent state, i.e., \(Sj)\ = 1. The values r\ £ {0, n/2} 
correspond to spin-nematic states with (Sj) = 0. For 
r/ = 0, all directors point along the z-axis (ferro-nematic 
state); whereas, for 77 = 7r/2, the directors on nearest- 
neighbor sites lie in a common plane at an angle of 120° 
(120° nematic state). 

As a second ordering pattern, we consider spin-nematic 
(NEM) states with (ip p \Sj -,\ip p ) = 0. The angle be- 
tween the directors on different sublattices is constant 
and taken as a variational parameter ( "umbrella" config- 
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FIG. 1: Variational energies (per site) for the bilinear- 
biquadratic model, Eq. (114[) . as a function of K, for D — —0.4. 
The system is iV = 12 x 12 lattice sites. 

uration). More precisely, we take the following family of 
spin-nematic states, 

d jeA = (0,- sin??, cos 77), 

, V3 . 1 . . (16) 

djeB.c = \^r— sm7 7, 2 sin 77, cos 77) , 

where the variational parameter 77 controls the angle be- 
tween the nematic directors on different sublattices. As 
before, the special value 77 = corresponds to a ferro- 
nematic, while 77 = tt/2 is a 120° nematic state. At the in- 
termediate value sin 77 = \/2/3, the directors are perpen- 
dicular to each other on neighboring sites (90° antiferro- 
nematic stated). 

Results 

Our variational results confirm the known phase dia- 
gram at D = 0. Furthermore, we find that the three- 
sublattice ordering of the ground state persists for non- 
zero values of D, i.e., all QSL states are higher in energy 
than the three-sublattice ordered states we considered, 3 - 
A typical plot of the variational energies (for D = —0.4) 
is shown in Fig.[TJ For K < 0.3, the magnetic Huse-Elscr 
wave function (J"- AFM), Eq. (JTUJ), is the best variational 
state. For 0.3 < K < 1, the fermionic antiferromag- 
netic state (f-AFM), Eq. ©, is the state with the lowest 
energy. As discussed above, the corresponding product 
states [specified in Eq. (|15l) ] are magnetically ordered 
with partially developed spins at 120° angles between 
sublattices. For D > 0, the ordered spins lie in the x-y 
plane while for D < 0, the spins order in a plane that con- 
tains the z-axis. For K > 1, the fermionic nematic states 
(f-NEM) take over, with directors specified in Eq. (|16l) . 
For D = 0, the best state is the 90° antiferro-nematic 
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stated and for D 7^ 0, the three nematic directors close 
(D > 0) or open up (D < 0) relative to the z-axis, de- 
pending on the sign of the single-ion anisotropy. In the 
fermionic long-range ordered states © , the optimal vari- 
ational parameter is h ~ 1.5. For this parameter value, 
the spinon excitations are fully gapped. Therefore, the 
spinons are confined and we expect that spin-wave ex- 
citations for these ordered states capture the low-energy 
physics of this mode L 12 ' 29 

The best quantum spin liquid states are the p+ip-state 
for K < — 1 and D ~ 0, and the unpaired U(l) state 
for K ~ 1, both having zero flux through the triangles 
(s = —1). All the other QSL states show very small or no 
condensation energies with respect to the £7(1) state. It 
is remarkable that, for K ~ 1, the U(l) state with three 
spinon Fermi surfaces is actually lower in energy than the 
optimized Huse-Elser wave function. We do not find any 
pairing instability of the U{\) spin liquid on the line K ~ 
1 for D > —0.8. For D < —0.8, there is a small energy 
gain from pairing in the d+id channel. However, the 
ordered fermionic states are still lower in energy. When 
D = 0, the three spinon Fermi surfaces match exactly. 
For D > 0, the f z Fermi surface expands while f x and f y 
Fermi surfaces shrink. The opposite happens for negative 
D. The kink in the U(l) energy in Fig. [1] marks the 
polarization to a ferro-nematic state with (S%j) = 0, for 
K < —0.6. That is, the spinon Fermi surfaces disappear 
at this point. 

The variational energies for the spin-one Heisenberg 
antifcrromagnet (K = D — 0) are displayed in Table fl] 
Note that the Heisenberg energy for the optimal product 
state of fully developed (coherent) spins ordered at 120° 
is —1.5. At the Heisenberg point, the spin liquids are even 
higher in energy than this uncorrelated product state. 

At the point K = 1 and D = 0, the model (TH]) enjoys 
an SU(3) symmetry^ 9 - (see also next subsection). On the 
line D — and arbitrary K, the remaining symmetry is 
SO(3) spin-rotation. This symmetry is broken to U(l) 
(generated by S z ) when D 7^ 0. However, on the line 
K = 1 and arbitrary D, the model possesses another 
SU(2) symmetry generated by the operators S z , S x — S* 2 , 
and S x S y + S y S x . The generator S x — S 2 allows to rotate 
the antiferromagnetic and the nematic ordered states into 
each other and they are degenerate. This property of 
the product states remains valid after the introduction 
of quantum fluctuations via §§§ or (|10[) . and it explains 
the degenerate crossings for the ordered states seen in 
Fig. [T] at K — 1. See Appendix iDl for a more detailed 
discussion of these symmetries. 



B. SU(3) ring-exchange model 

In the last subsection we concluded that the simplest 
extension of the spin-one Heisenberg model, Eq. (fl4|). 
does not show quantum spin liquid behavior. To moti- 
vate another promising spin-one model, let us consider 
an SU(3) symmetric Hubbard model for three flavors of 



Variational state 


Heisenberg energy 


Huse-Elser J-AFM 


-1.783(1) 


Fermionic f-AFM 


-1.570(2) 


p+ip spin liquid 


-1.33(0) 


(7(1) spin liquid 


-1.00(3) 



TABLE I: Variational energies (per site) for the spin-one 
triangular-lattice Heisenberg antiferromagnet, (I14|l . for K — 
D = 0; N = 144 sites. 

fermions / a , 

Hsu(s) = "* E + U H n h ( 17 ) 

where nj = J2 a n a,j = J2 a fajfaj- Let us consider the 
case when each flavor is at 1/3-fhling (J^j n a j/N = 1/3). 
For U 3> \t\, the low-energy subspace of this model 
corresponds to the spin-one Hilbert space. Similar to 
Refs. (mm, we can derive a low-energy effective spin- 
one Hamiltonian for (|17p . To lowest order in t, we find 
the exchange term 

E faifbifbjfaj ■ (18) 

To next order, the following three-site term is expected 
to arise: 

E {flifbifUvflfak + h.C.} , (19) 

where the sum k) is over elementary triangles of the 
lattice. Let us write the flavor exchange operators in (fT5)) 
as Vij = ^abfaifbifbjfaj- The three-site terms in ([19]) 
correspond to PijPjk +'Pjk'Pij ■ These operators move the 
local states clock- and anticlockwise around the triangles 
of the lattice. 

In the case of a similar Hubbard model with two 
fcrmion flavors (spin S = 1/2), the exchange operator 
Vij appearing in the low-energy model corresponds to the 
Heisenberg term in spin language, Vij — 2Si ■ Sj + 1/2. In 
this case, a three-site term VijVjk is trivial in the sense 
that it can be reduced to a sum of two-site terms. For 
three flavors, however, the situation is different. In that 
case and spin S = 1, one finds^ 

Vij - S, ■ Sj + (Si ■ S^ 2 - 1 . (20) 

Therefore, the lowest-order term (fT5|) corresponds to the 
KD-model (HD with K = 1 and D = 0. The next- 
order ring-exchange term (fT9|) is a non-trivial perturba- 
tion since it cannot be reduced to two-site terms. Ring- 
exchange models for spin 1/2 on the triangular lattice 
with non-trivial four-site plaquette terms^i are believed 
to exhibit spin-liquid ground states i 38 ' 39 

Motivated by the three-flavor Hubbard model, we pro- 
pose to study the SU(3) symmetric ring-exchange model, 

H a = cos a *Pij + sm a E {'Pn'Pik + h - c -} ■ ( 21 ) 

(*ii> (h3,k) 
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FIG. 2: Pictorial presentation of the variational phase dia- 
gram that we find for the SU(3) ring-exchange model (|21[) . 



■U(l) QSL 

■ d+id QSL 

■ JT-NEM 

■ f-AFM 
FM 




FIG. 3: Variational energies (per site) of the SU(3) ring- 
exchange model, Eq. (|21|l . as a function of a/n. N = 12 X 12 
lattice sites. 



The sum in (|21|) goes over nearest-neighbor links 
and elementary triangles (i,j,k) of the lattice. The pa- 
rameter of this model is a € [— 7r,7r]. 



Results 

An analysis of the model (f^Tj) in terms of general three- 
sublattice product states reveals a ferromagnetic phase 
in the parameter range a < — arctan(3/4) ~ — 0.2-7T 
and a > n — arctan(l/2) ~ 0.857T. Any uniform prod- 
uct state Y[j \ a )j i s an exact eigenstate with energy 
e(a) = 3 cos a + 4 sin a, and it is the lowest-energy three- 
sublattice product state in this parameter range. For 
- arctan(3/4) < a < arctan(3/2) ~ 0.31tt the 120° anti- 
ferromagnetic product state is stabilized. Finally, in the 
range arctan(3/2) < a < it — arctan(l/2), the three ne- 
matic directors order in a common plane at an angle of 
120° to each other on nearest-neighbor sites, (see Fig. [2]) 

Above analysis with uncorrelated three-sublattice 
product states revealed the same ordering patterns as 
we found in the case of the K D-mode\ investigated in 
the previous subsection. Therefore, we can use the same 
trial wave functions specified in Eqs. (|15l) and (|TB)) to 
construct correlated ordered states for the ring-exchange 
model. 
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FIG. 4: Optimized variational parameters A = \A xy \ (dot 
symbols, left scale) and (S'f) — 2/3 (x symbols, right scale) for 
the d+id QSL state in the ring-exchange model (|2ip . Among 
the states we consider, the d+id state has the lowest energy 
in the range 0.17tt < a < 0.33tt. For 0.17tt < a < 0.22tt, the 
optimal state is a 0-flux state with s — — 1; for 0.22-7T < a < 
0.337T, we find a 7r-flux state with A ~ 0.5 and s = 1. 



We calculate the variational energies of the QSL states 
([7]) as well as the energies of correlated three-sublattice 
ordered states ([9]) and (jTO)) for the ring-exchange model 
specified in Eq. (|2"T|) . The results are presented in Fig. [3] 
We see that the conclusions we draw from the simple 
product state calculation above agrees with the result 
using correlated wave functions in most of the param- 
eter range. However, in the region between the AFM 
and the 120°-nematic phase, around a ~ 7r/4, we find an 
extended region where the d+id QSL has the lowest en- 
ergy (see also Fig. [2] for a scheme of the phase diagram). 
The optimal d+id variational parameter lA^I along with 
(5^> -2/3=1/3- N z /N are shown in Fig. GO The ring- 
exchange term favors a 7r-fiux d+id state with s = 1: As 
a increases, the 0-flux state with a large pairing term 
(| A*" | ~ 4) changes to a tt-Aux state with |A^| ~ 0.5 at 
a ~ 0.22tt. 

Note that the d+id QSL phase in Fig. [2 as well as the 
adjacent 120° nematic phase, exhibit a ferro-quadrupolar 
order. For the d+id state this is apparent from Fig. [H 
since (S^) > 2/3. In contrast, lattice rotation symmetry 
is unbroken in the d+id QSL while both adjacent ordered 
phases spontaneously break lattice rotation. 

As discussed above, for a = and up to a constant, 
(j2"Tj) corresponds to the model (jT4j) with K — 1 and 
D = 0. The ground state of this model was recently 
approached with density matrix renormalization group 
(DMRG) calculations in Ref. [ioj]. In this work, the au- 
thors found a three-sublattice ordering pattern that is 
consistent with our result. The DMRG energy is dis- 
played in Table |TT] along with the variational energies of 
the lowest-energy states used in the present paper. 
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State 


5(7(3) energy 


Fermionic f-AFM 


-0.57(8) 


Huse-Elser J7-AFM 


-0.27(7) 


(7(1) spin liquid 


-0.34(3) 


DMRG-^ (jV = 8 x 10) 


-0.678 



TABLE II: Variational energies for the 5(7(3) model, Eq. J21 
at q = on JV = 12 x 12 sites. 



We also consider additional perturbations to the ring- 
exchange model (|2"Tj) in order to assess the effect of such 
terms on possible low-temperature QSL phases. First, 
we add a single- ion anisotropy term D J2j ■ As dis- 
cussed in the previous section, such a term breaks the 
5(7(3) symmetry of the model to 5(7(2). For small D, 
the ordering plane is explicitly chosen. Large D deforms 
the three-sublattice ordering pattern in a similar way as 
in the bilinear-biquadratic Heisenberg model (jl~4)) . How- 
ever, we find that the phase boundaries of the d+id state 
with the adjacent ordered states are barely affected by 
D. (We investigate the range \D\ < 1.5) 

Second, we add a next-neighbor exchange term 
^E{{ifl)^3 to (p?Tj). Such a term strongly frustrates 
three-sublattice ordering. At a = 0, we find that the 
three-sublattice ordering is destroyed for J2 as small as 
J2 ~ 0.25, and the (7(1) QSL has the lowest energy 
among our ansatz wave functions. However, an analy- 
sis of this model in terms of product states reveals that 
the competing ordering pattern is spiral. So far, we have 
not included spiral states into our variational analysis 
and we reserve a detailed study for future work. 

In conclusion, we find a region in the phase diagram 
of the 5(7(3) ring-exchange model (|2"Tj) where a chiral 
d+id quantum spin liquid state is stabilized. The ques- 
tion remains whether this spin model can describe the 
relevant magnetic interactions in the 6H-B structure of 
Ba3NiSb20g. A perturbative expansion in t/U of a two- 
band Hubbard model with an additional orbital degree 
of freedom and strong Hund coupling would produce a 
spin 5=1 Heisenberg term Si ■ Sj to order t 2 . Only to 



next order, i 4 , one expects bi-quadratic terms (Si ■ Sj) 2 
as well as three-site terms (Si ■ Sj)(Sj ■ Sk)~~— Those 
three-site terms are present in our ring-exchange model, 
but we need them to be of the same order as the Heisen- 
berg term. We found that a dominant nearest-neighbor 
Heisenberg term is detrimental to the stability of the 
d+id quantum spin liquid. Further terms in (|2"Tj) like 
(Si ■ Sj) 2 (Sj ■ Sk) would only come to order t 6 or higher 
in a perturbative expansion. At present, we cannot con- 
clude if these higher order terms are relevant for the sta- 
bility of the d+id QSL or not. 

While it is unclear, at present, if the ring-exchange 
model (f2"Tj) can realistically describe the spin-liquid phase 
of Ba3NiSb2 0g, it is a very natural model to study 
if one starts from an integer-filled three-band Hubbard 
model. Such three- (and higher-) band Hubbard mod- 
els are currently of great interest, both theoretically and 



experimentally, in the cold atoms community; see, e.g., 
Refs. EUl. 



V. CHIRAL EDGE MODES FOR THE d+id QSL 

In Ref. (45|, it was shown that the d+id superconduc- 
tor is a topological state with Chern number equal to 
two. From the bulk-edge correspondence, this indicates 
the presence of two chiral edge modes. A semiclassical 
argument^ supports this conclusion. Let us first recapit- 
ulate the semiclassical argument and generalize it to chi- 
ral topological superconductors. Next, we discuss these 
results in the framework of the d+id QSL state. 



A. Edge modes in chiral superconductors 

In the bulk of a superconductor involving two fermion 
flavors, writing tp = (ci, c|), the Bogolubov equations are 




i'k = . 



(22) 



Here, we consider fully gapped superconductors with 
|A(fc)| > 0. The spectrum is given by Ek = 
±^ 2 + \A(k)\ 2 . 

Next, consider a superconductor with a boundary 
along the x-direction. Asymptotically (i.e., for \k\y 3> 1), 
an incident bulk wave packet ipk,e lk r with k = (k x ,k y ) 
is reflected at the boundary to an outgoing wave packet 
ipkie lk r with wave vector k 1 — (k x , —k v ). The two wave 
packets "see" the gap functions A(fc) and A(fe'), respec- 
tively. For a given incident wave vector k, it seems there- 
fore possible to map this problem on the half plane to a 
one-dimensional scattering problem where the order pa- 
rameter A(y) interpolates from A(fe) as y —¥ —00 to 
A(fe') as y — >• +00: 



id y - E A(y) 
id v — E 



A(y) 



(23) 



This one-dimensional problem can now be solved in the 
usual way 47 ' 48 For simplicity, let us choose a bulk po- 
tential of the form A(fc) = Ae m ^> where I e Z is 
the winding of the phase of the order parameter around 
the Fermi surface, and cos6(k) — k x /\k\. A scattering 
state with incident angle 9 results in an outgoing angle 
9' = ir — 9. The asymptotic potentials can therefore be 
chosen as A(fc) = Ae"* 9- */ 2 ) and A(fc') = Ae-^ 6 ^^. 
Accordingly, the order parameter A(y) in (|23[) has a 
constant real part Acosl(9 — tv/2), and an imaginary 
part A sin 1(9 — n/2) that changes sign across the bound- 
ary. This problem can be solved exactly for certain spe- 
cial cases of the interpolating gap profiles ! 48 ' 50 For ex- 
ample, a kink profile with A(y) = A[cosZ(0 — tt/2) — 
i tanh(y) sin 1(9 — tt/2)] yields the bound-state spectrum 



E e = Acos[l(9 - -)] sgn[sinZ(0 - -)] 



(24) 
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FIG. 5: The four lowest energy levels of the d+id mean-field 
state (|7|) on an infinite triangular-lattice strip as a function 
of wave vector k x along the strip. The width of the strip is 
200 sites. The boundaries are chosen to be parallel to one 
lattice direction and we use open boundary conditions. The 
spectrum of the f z spinon with a bulk Fermi surface is omit- 
ted. The gapless states (blue online) are localized on the 
lower boundary for left movers (dashed line), and on the up- 
per boundary for right movers. The higher states (red online) 
are delocalized and the energy levels above them are "dense" . 



with corresponding eigenvector 

Mv) = w ; (l,sgn[sini(0- £)]) . (25) 
2 cosh(y) 2 

We observe that Eg as a function of 8 vanishes exactly \l\ 
times. Therefore, there are \l\ gapless edge modes. Using 
cos0 ~ k x /kf and expanding Eq. (pM)) around a node at 

momentum fc™, we get Eg ~ — l(k x — k x )A/kf + The 

edge modes are chiral and propagate with the velocity 
v n = —lA/kf. In the simplest nontrivial and well-known 
case of a p+ip-superconductor— ~— (Z = 1), a single chiral 
edge mode is located at fc™ = 0. Higher angular momenta 
have chiral modes at fc™ + 1 0. Since the phase winding 
of the order parameter around the Fermi surface is a 
topological property, we expect that the number of chiral 
edge modes is a robust feature of the state, too. On the 
other hand, the precise location of the nodes {A:™} and the 
corresponding propagation speeds \v n \ depend on further 
microscopic details. 



B. Low energy edge theory for the d+id QSL state 

According to the above semiclassical argument, the 
d+id QSL state (I = 2) is expected to exhibit two chi- 
ral edge modes located at wave vectors k x ~ ±fc//\/2. 
To substantiate this claim, we calculate the spectrum of 



the d+id state (JT)) on an triangular-lattice strip of in- 
finite length [Here, we neglect the local constraint ^ 
and work in the fermionic Fock space]. The four lowest 
energy levels are shown in Fig. [5] as a function of wave 
vector k x along the strip. The triangular-lattice d+id 
state indeed exhibits two gapless left movers localized on 
the lower boundary and two right movers localized on 
the upper boundary. The spectrum of f z spinons with a 
bulk Fermi surface is omitted in Fig. [5] 

As discussed above, the low-energy degrees of freedom 
localized on the edge for the d+id QSL state are two 
chiral Dirac fermions. To discuss the physics of these 
edge states, it is convenient to go to the spinon basis 
creating S z eigenstates. We have 

fa = ^(fx-i(Tf y ), (26) 

with a — ±1. We also denote /j = /_i . The x-y (triplet) 
pairing term of the d+id state is fxifyj — fyifxj — 
i(fnflj ~ fnfij)- We consider an edge along the re- 
direction and denote the momentum along the edge by 
k = k x G [— 7t, 7r]. The two gapless points in the bound- 
ary spectrum are denoted by fc™ = ±fc° with k° > 0. 

Using the semiclassical expression (|25"|) , the edge states 
are created by operators 

X*(k)~ f ak +asgn(k)fl_ k , (27) 

for \k\ ~ k°. The excitations Xi(fc) and Xl(k) carry spin 
S z = ±1, respectively. Note that the edge states at pos- 
itive and negative momenta k are not independent: We 
have Xa(k) = ersgn(fc) \s The low-energy effective 
edge Hamiltonian is given by 

n = v (k-k°)xUk)xAk), (28) 

where the sum over k is restricted to the vicinity of 
the node at momentum +fco to avoid double counting 
of states. 

Similar to the ordinary quantum Hall effect, the chi- 
ral edge modes are expected to be robust with respect 
to disorder and impurities because no backscattering 
is possible^ Furthermore, due to S z conservation, hy- 
bridization terms like f\x<y cannot appear in the low- 
energy Hamiltonian. The presence of chiral edge modes 
carrying spin S z — ±1 implies a quantized spin Hall con- 
ductivity. We also expect a thermal Hall conductivity in 
the d+id QSL state. 

In the d+id QSL phase with unbroken lattice symme- 
tries, f z must necessarily form a spinon Fermi surface. 
However, this argument becomes invalid when the lat- 
tice symmetries are explicitly broken. For example, close 
to the edge of the sample, the f z spinon is susceptible 
to pairing. Similarly, we expect the spinon to acquire a 
local gap in the vicinity of bulk impurities. This prop- 
erty makes the spinon Fermi surface hard to detect in 
any experiment that involves local probes. 
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VI. RESPONSE FUNCTIONS AND PHYSICAL 
PROPERTIES OF THE d+id QSL 

In this section, we release the local constraint ([3]) in 
order to analyze the spectral properties of the d+id mean 
field state. This can be justified from the point of view of 
the U(X) gauge theory since we are in a Higgs phase where 
gauge fluctuations can be neglected. In this case, the 
f z spinon can be treated as a weakly interacting Fermi 
liquid. 



Static spin susceptibility and NMR relaxation 
rate 



The 



response 



function, 



Raa{iu) 



Eij/o dr e lUT {S ai {T)S a j), in the d+id state has the 
following properties: Since f x and f y fermions are paired, 
we have R zz {iuj) — 0. However, R xx (iuj) = R yy (iuj) do 



not vanish at low temperature, 
low-temperature limit, \uj\ <C T - 



In the low-frequency, 
• 0, we find 



xl = Re[R xx (0)} = 



cPk E k - sgn(ff )gg 
BZ 2tt Ek (E k + \a\) 



(29) 
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FIG. 6: Wilson ratio, (|31|) . for the d+id state as a function 
of the spinon chemical potential, fi z — fi%. The shift cor- 
responds to the optimal value of the chemical potential in 
the ring exchange model (|21[) at a = 7r/4 (without single-ion 
anisotropy) . 



where £g = 2s[cos(l • k) + cos(2 ■ fe) + cos(3 • ft)] 
is the dispersion of the x- and z-fermions. E k = 
V (£fc) 2 + l^fe W | 2 5 an d the d+id gap function is A^ = 
A[cos(i-ft)+e 27ri / 3 cos(2-fc)+e- 27ri / 3 cos(3-ft)]. As before, 
1, 2, and 3 are vectors of nearest-neighbor links on the tri- 
angular lattice. We find that the static spin susceptibility 
X x takes a non-zero value given by the integral over the 



Brillouin zone (BZ), Eq. (|29|) . Its numerical value de- 
pends on the parameters A, /j, x , fi zi and on s = ±1. 



In the limit A < 



A. x 



approaches the 



Pauli susceptibility of two unpaired Fermions, x x — 2^ z , 
where v z = f BZ d 2 k/(2n) <5(£|) is the density of states 
at the Fermi surface. In conclusion, we predict a strong 
anisotropy of the spin susceptibility xl = Re[i? O a(0)] for 
the d+id state at low temperature. 

The nuclear spin relaxation rate is given by T^ 1 ~ 
Tlm[R(iuj -> 0)]. In the d+id QSL state, we find that 
this quantity is exponentially small for temperatures be- 
low the gap. 



B. Specific heat and Wilson ratio 

In the d+id spin liquid, the magnetic specific heat at 
low temperature is linear in temperature due to the f z 
spinon Fermi surface. The coefficient is given by£2 
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The Wilson ratio is defined as 

4tt 2 xo 
3 7 



Rw 



(30) 



(31) 



Since the measurements in Ba3NiSb20g were made on 
powder samples, a directional average should be used in 
this expression for comparison with experiment, xo = 
2*2/3. 

The Wilson ratio, Rw = 8%2/(3^z), for the d+id state 
is plotted as a function of \x z — n° z in Fig. [6] The choices 
of parameters (A = 0.5 and 2.6 for the 0-flux state, and 
A = 0.5 for the 7r-flux state) are examples of lowest en- 
ergy d+id states in the ring-exchange model, (|2"T1) , at 
a ~ 7r/4. Note that, in this plot, we adjust the chemical 
potential [i x such that the constraint is satisfied on aver- 
a 9 e i 2 a ( n<1 ) = 1- The shift [i z is the optimized chemical 
potential for the ring-exchange model, i.e., for D = 0. In 
Fig. [6] we see that the Wilson ratio is enhanced in the 
d+id state with respect to a metal (where Rw = 4/3) 
by a factor of approximately two at fi z = /j,®. This can 
be attributed to the fact that only a single fermion flavor 
contributes to the coefficient of specific heat in the QSL 
state. Since we have D oc {jjl z — /i^), Rw can be further 
enhanced by a single-ion anisotropy term in the Hamil- 
tonian. An easy-plane anisotropy (D < 0) shrinks the 
spinon Fermi surface, resulting in enhancement of Rw- 
For an easy-axis anisotropy (D > 0), the Wilson ratio is 
enhanced due to an increase in magnetic susceptibility in 
the case of the 7r-flux state. 

Experimentally, a large Wilson ratio of Rw — 5.6 
was reported for the spin liquid phase of Ba3NiSb20g. 11 
Within the framework of the d+id QSL state, we can con- 
clude that quite a strong single-ion anisotropy, \D\ > 1, 
is required to explain the large Wilson ratio seen in 
BagNiSbsOg. 

Note that the U(l) state has Fermi surfaces for all 
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three spinons. However, since this state is in a Coulomb 
phase, the f/(l) gauge fluctuations are expected to be 
very strong. Assuming Landau damping of the photon, 
it has been proposed that the specific heat in such a 
strongly coupled phase should have the non-Fermi-liquid 
behavior C M oc T 2 / 3 .^^ 



the Cu 2+ -ions on the triangular lattice may form dipo- 
lar molecules with the Sb + -ions and can move out of 
plane. Strong disorder due to Jahn- Teller distortions or 
fluctuations of these Ising dipoles may play a key role 
in the absence of ordering in the Cu-compound. Similar 
effects may also be present in Ba3NiSb2 0g, which opens 
promising avenues for future studies on this material. 



Thermal Hall effect 



Due to the spinon Fermi surface of f z , the d+id QSL 
state exhibits a longitudinal heat conductivity^ 3 - Accord- 
ing to the Wiedemann-Franz law, it is of the form 



9o ■ 



(32) 



where go = tt 2 T / (3K) is the thermal conductance quan- 
tum, ef is the Fermi energy of the f z spinon, and t is 
its lifetime. However, no longitudinal spin current will 
flow since the spin excitations are fully gapped in the 
bulk. Nevertheless, we expect a thermal (and spin) Hall 
conductivity due to the chiral edge modes; 54 ' 55 



2<7o- 



(33) 



Since the state is compressible, n xy is not expected to be 
exactly quantized. The f z spinon with a bulk Fermi sur- 
face also contributes to K xy due to a classical Hall effect 
in the chiral spin liquid. On the other hand, the spin 
Hall conductivity is expected to be exactly quantized. 
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Appendix A: Variational Monte Carlo 

The Variational Monte Carlo (VMC) method allows 
to efficiently evaluate expectation values of observables 
in a given many body wave function within small error 
bars . 22 ' 57 This works as follows: Let be the wave 
function and let O be the observable we want to evaluate. 
Let {\a}} be an "Ising" basis of the Hilbert space, i.e., 
| a) is product of local basis states. We can write 



\o\rp) 



(Al) 



VII. CONCLUSION AND OUTLOOK 

In this paper, we construct all natural quantum spin 
liquid states with three flavors of fermionic spinons for 
spin S = 1 Heisenberg models on the triangular lattice. 
We compare their variational energies with the ones of 
various long-range ordered states. We find that, for large 
biquadratic and ring-exchange terms (of the order of the 
Heisenberg exchange J > 0), an exotic chiral quantum 
spin liquid with a spinon Fermi surface is stabilized. The 
physical properties of the d+id QSL state seem to be 
consistent with the recent experiment on Ba3NiSb2 0g. 10 

While the d+id QSL scenario we investigate in this pa- 
per has many attractive and novel features, it remains un- 
clear if the microscopic parameters required to stabilize 
such a phase are realized in Ba3NiSb2 0g. From the crys- 
tal structure proposed in (Tlj . it seems more likely that 
the nearest-neighbor antiferromagnctic exchange energy 
J is the dominant microscopic parameter. Therefore, the 
theoretical research must continue and more experiments 
are needed to elucidate the spin state realized in this ma- 
terial. 

Recently, new experimental results were published on 
the related spin- liquid candidate Ba 3 CuSb 2 0g in [56| . 
In contrast to earlier experiments on powder samples^ 
the new experiments on single crystals indicated that 



with |^(«)| 2 = |(o#)| 2 /(^). Since £ a \^(a)\ 2 = 1, 
| ip{ a ) 1 2 is a probability distribution on the Ising config- 
urations {a}. Such a distribution can be generated by a 
Metropolis algorithm with acceptance probability 



a = mm 



ip(a) 



(A2) 



Note that, in (|A2|) . -0(a) does not need to be normal- 
ized. The sequence {a} generated by a random walk 
with probability (|A2I) can be used to efficiently calculate 
the expectation value, 



(A3) 



In this paper, we use ~ 200 Monte Carlo runs to es- 
timate the error of Eq. (|A3[) by its variance over the 
runs. The length of a run is ~ 200 steps, and the observ- 
ables are measured after each step. The measurements 
are precessed by an equilibration skip of ~ 400 steps. 
Each Monte Carlo step consists of 4 x N ~ 600 local 
moves, accepted with probability (|A2I) . We use a lattice 
of N = L x L sites, with a linear size L = 12 in our 
calculations. The error bars on the variational energies 
shown in Figs. [T] and [3] are smaller than the symbol sizes. 
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Local and global constraints (projections) on the wave 
function \tp) can be easily implemented in the VMC 
scheme. The Gutzwiller projection, \ip) — Pq\iPo), for 
example, can be taken into account by restricting the 
Ising configurations a to the singly-occupied subspace 
(rij = 1). Similarly, projection of a spin wave function 
to 5*°* = leads to a global restriction on the configura- 
tions a. Here, it is important to have an algorithm that 
generates all states a in the constrained subspace with 
uniform probability. 

To apply VMC to a particular wave function, we first 
need an expression for tp( a ) ( a \4'}- Next, an effi- 
cient algorithm is needed to calculate the Metropolis ac- 
ceptance probabilities (|A2j) for local moves in the con- 
strained subspace. Similarly, for each observable of in- 
terest, one has to find an efficient way to calculate the 
ratio of overlaps in (|A3I) . 

Appendix B: Fermionic wave functions 

The first class of wave functions that we are con- 
sidering in this paper are Gutzwiller-projected ground 
states of quadratic Hamiltonians, Hmf, for three flavors 
of fcrmions f a . A similar study of wave functions with 
two flavors of fermions has been pioneered by Gros^ 2 - for 
spin S = 1/2 models. 

In our calculation of fermionic QSL and fermionic or- 
dered wave functions, we use the local basis of time- 
reversal invariant states, | a) € {|x), \y), \z)}, Eq. @. The 
Ising configurations a are restricted to singly occupied 
states on a lattice of A" = L x L sites. Furthermore, we 
restrict the configurations to states with N x = N y and 
A" 2 kept fixed (That is, the wave functions are projected 
to fixed total flavor numbers; see below). 

Let r° e Z L x 7L L be the lattice positions of flavor 
a G {x, y, z} in the Ising configuration a. The U(l) state 
and the triplet (x-y paired) QSL states in ([7]) can be 
written as a product of two determinants, 

(a\il>) = det[e ik r r i] detL4(r* - if)] , (Bl) 

where j and I are the indices for the determinants, fcj 
are the occupied momentum states of f z spinons inside 
the Fermi sea, < fi z . For the U(l) state, A(r) is a 
Slater matrix^ 

A(r)= elk r > ( B2 ) 

fceBZ, 

with momenta k going over filled states in the first Bril- 
louin zone (BZ). For the triplet QSL states (s-wave, 
d+id), we have 

Mr)= ]T a k e lkr , (B3) 
fceBZ 

where a k = Vk/uk — Ak/(Ek + £&) is the ratio of BCS 
coherence factors for the pairing of f x and f y fermions. 4 



For the QSL states with equal-flavor pairing (/- 
wave, p+ip), the wave function is a product of three 
Pfaffians, 5 ^ 9 . 

(a\^) = Y[Pi{A a (r«-r?)}, (B4) 

a 

with 

A a (r)= J2 a k sin(fc-r), (B5) 
fceBZ 

where = v^/u^ are the ratio of coherence factors for 
each paired fcrmion flavor. 

In the case of the ordered states ([9]), the fermions are 
unpaired, but the flavors hybridize through terms f} a fjb, 
etc. For a lattice of N sites, the corresponding wave 
function is a single Slater determinant of size N x N, 

(a\1>) = det[A i (r°)]- (B6) 

Here, I = 1...N, and Ai(r°") are the lowest eigen- 
vectors of the mean-field matrix Hf^ with i? or d = 

ab fai^ijfbj [F° r the three-sublattice ordered states 
we consider in this paper, the eigenvectors can be labeled 
by I = (n, k), where n is a band index and k lies in the 
reduced Brillouin zone]. 

Our calculations are done on a finite lattice with A^ = 
Lx L sites. In order to avoid singularities or degeneracies 
in (|Bip - (|B6[) . we use quadratic trial Hamiltonians ([7]) and 
§§§ with periodic in one, and antiperiodic boundary con- 
ditions in the other lattice direction for the spinons f a j. 
The /-wave state, however, has lines of nodes in the gap 
function (at momenta {fco}) that cannot be avoided 
by choosing periodic-antiperiodic boundary conditions. 
A singularity |a£ | — > oo occurs on these lines, and (|B5|) 
is ill-defined. To cure the divergencies, we replace a^ o by 
a large but finite quantity, namely, ±20 x max fc ^{ fco }|a^|. 
The sign is chosen to be consistent with the sign of a£ as 
k — > feo- We have verified that the relevant correlators do 
not depend on the precise factor in the regularization and 
that the wave function (correlators) correctly reproduces 
the U(l) state when |A aa | < 1. 

We use the usual tricks for an efficient evaluation of the 
Metropolis acceptance probability (|A2[) and the expecta- 
tion values (|A3p in fermionic wave functions: The inverse 
of the matrices in (|F31[) , (|B4|) , and (|B6|) is stored and up- 
dated during the Monte Carlo random walki^i This al- 
lows for efficient evaluation of determinants and Pfaffians 
with rows and/or columns replaced or removed. 22 ' 59 ' 60 To 
update the inverse of an anti-symmetric matrix with a 
row and column replaced, we use the Sherman-Morrison 
algorithm 6 ^ twice, followed by anti-symmetrization of the 
matrix. This procedure greatly improves the numeri- 
cal stability of the update. The "pfapack" package by 
Wimmer— is used for efficient evaluation of Pfaffians. 
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Flavor-number nonconservation 

An important technical difficulty with fermionic RVB 
wave functions for spin S = 1 is that typical microscopic 
models like (1141) . when written in terms of fermion op- 
erators, do not conserve the number of each fermion fla- 
vor separately. This issue is also present if we wish to 
represent the spin operator by more than three fermion 
flavors. We have N a = N - J2j S 2 aj and [H KD ,N a ] ± 0, 
in general. Unlike in the case of spin-1/2, conservation 
of 5* ot = J2j Szj does not imply conservations of fla- 
vor number. Note, however, that N a is conserved in 
the SU(3) model, (HQ, or in the KD-mode\ (HU at 
K = 1, where this issue does not arise. Writing (fT4")l 
with fermions, the terms not commuting with N a are 



(B7) 



ab 



which vanish for K = 1. In general, there is therefore 
no justification for using variational wave functions that 
are particle-number eigenstates. For such wave functions, 
the Ising configurations a in (IA2|) must visit all possible 
total flavor numbers, with ^ a N a = N kept fixed. In a 
brute force implementation, the determinants and Pfaf- 
fians in (|B1[) and (|B4[) may need to change sizes during 
a Monte Carlo run, which implies a high computational 
overhead. Such a simulation has recently been done in 
the case of spin-one chains^ 

The problem is actually absent for the QSL states with 
a spinon Fermi surface. In this case, the wave function 
is an N z eigenstate. N x and N y do fluctuate in a paired 
state; nevertheless, N x = N y and all expectation values 
of (IB7|) vanish in this class of wave functions. The dif- 
ficulty is only present for equal-flavor paired QSL states 
(/-wave and p+ip) and for the ordered states ©. In 
these cases, the expectation value of (|B7|) docs not vanish 
(before or after Gutzwiller projection). The flavor num- 
bers N a fluctuate independently of each other in these 
wave functions. 

To resolve this issue, we can use the standard 
argument^ that relates grand-canonical and micro- 
canonical RVB wave functions: The paired mean-field 
states are strongly peaked at some average flavor num- 
ber N = (N°,N°,N°). This peak in flavor number 
may shift position to No after Gutzwiller projection, 
but it should still be present. Furthermore, the vari- 
ance is expected to vanish in the thermodynamic limit, 
((N a - N°) 2 )/N 2 ~ 1/JV. Therefore, it is justified to 
work with micro-canonical wave functions that are ob- 
tained by projecting the grand-canonical wave function, 
\ip), to fixed total flavor numbers, 



\N ) = P(N )\if>). 



(B8) 



VMC calculation of expectation values of particle- 
number conserving operators within a micro-canonical 
wave function is straightforward. However, off-diagonal 



operators like (|B7[) require some care*^ As an example, 
let us consider the operator 

Rxy = fxifyifxjfyj ■ (B9) 

Its expectation value in the grand-canonical wave func- 
tion can be approximated as 



(N+\R xy \N Q ) 



(N+\N+)(No\N ) 



(BIO) 



with = (N° ±2,N°T 2, N°), and N is the average 
particle number in \ip). In VMC, the right-hand side of 
Eq. (|B10[) cannot be calculated directly with the correct 
normalization. However, it is possible to calculate 

(N+\R xy \No) ^ (JVqI^IJVq-) 



(N \No 



(N \N ) 



within a single Monte Carlo run. Since the 
last average satisfies |(JVo|i2a;y|JVJ")/(JVo|JVo) — 
| (Nq \R xy \N a ) I (Nq I JVq") J , the normalization factor can 
be calculated from the ratio of the two correlators 
in (Em, 



9xy — 



(Nq\Nq) 



(N \R xy \N - 



(N+\R xy \No) 



(B12) 



Finally, the correctly normalized expectation value (|B10j) 
is evaluated as 



(B13) 




(AWo) 



It is clear that, for a given wave function, g xy , (|B12I) . 
does not depend on the off-diagonal operator R xy (For 
example, R xy on different sites must give the same g xy ). 
This provides a nontrivial check of our code and we found 
that the renormalization factors g a b are indeed identical 
on different sites within error bars. 

Of course, a particle-number projection \N) is only 
a faithful representation of \ip) if the flavor number 
N is sufficiently close to the average value iVo in the 
Gutzwiller projected wave function. Using TV as a vari- 
ational parameter (here with the restriction N x = N y ) 
guarantees that the state |iVo) cx \ip) is among the varia- 
tional wave functions. For the equal-flavor paired singlet 
wave functions, we found that the agreement between our 
optimal correlators and the ones calculated in the corre- 
sponding grand-canonical wave functions is very good^ 

For spin S = 1/2 systems, the investigation of (doped) 
RVB wave functions in the grand-canonical ensemble was 
pioneered by Yokohama and Shiba in Ref. [66] . These au- 
thors introduced a particle-hole transformation c| i— > 
that allows to do fixed-particle VMC simulations. How- 
ever, this trick does not easily generalize to spin-one. 
For spin-half RVB wave functions, the agreement be- 
tween micro-canonical and grand-canonical approaches 
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was found to be very good. Note, however, that parti- 
cle number renormalization by the Gutzwiller projector 
in grand-canonical wave functions leads to subtle effects 
that need to be taken into account if one wishes to apply 
the Gutzwiller approximation^ - — 

Appendix C: Huse-Elser wave functions 

This appendix contains details regarding the imple- 
mentation of trial wave functions of Huse-Elser type, gen- 
eralized to the spin S — 1 case. Similar to the case of 
spin S = 1/2j22i our construction starts from an uncor- 
related product-state wave function. Quantum correla- 
tions are introduced by applying Jastrow factors to the 
simple product state. The resulting wave function has 
two sets of variational parameters: parameters control- 
ling the product-state, and Jastrow parameters respon- 
sible for the quantum correlations. 

For the Huse-Elser wave functions, we use the local ba- 
sis of S z eigenstates, i.e., the states |0), |l), and |l) with 
S z = 0,1, and — l, respectively. The corresponding basis 
of Ising configurations is denoted by \a) = 1 1 101011 .. .}. 
As before, the singly occupied subspace corresponds to 
physical spin states. Furthermore, we project the wave 
functions to iS* ot = by restricting to Ising states with 
Ni = Nj. However, here we allow the total flavor num- 
bers to fluctuate within this subspace. 

In Ref. [3l| the optimal three-sublattice product states 
for the bilinear-biquadratic model (|T4|) were calculated. 
It was found that the ordering patterns in this model 
are well captured by the antiferromagnetic and nematic 
states given in Eqs. (TT51) and (p~6|) . In the basis of S z 
eigenstates, the wave function on A, B, and C sublattices 
is given by 

|A)=cosr ? |0>+ K ^(|l) + |l)) ) 

\B), \C) = cost, |0) + k-—^=J- (e^\l) + e ±2 ^ |1» , 



where 77 is a variational parameter, k, = 1 corresponds 
to the antiferromagnetic, and k = i to the nematic prod- 
uct state. Using the Ising basis, the corresponding wave 
function may be written as 

IW}=£e"», (C2) 

a 

where the sum goes over Ising states in the S z basis. 
The one-body operator Hi accounts for different weights 
of |0), |1), and |1), as well as for site-dependent phase 
factors in the product-state wave function. For the par- 
ticular case specified in Eq. (|Clj) . it can be written in 
terms of S z operators as 

Hi =Y J {^^c^ ieB )S Z] +\og{n-^-)S%}. (C3) 



The Kronecker symbols 5j^B and Sj^c are non-zero only 
for sites j belonging to the B or C sublattice, respectively. 

The advantage of the rather complicated form (|C2j) for 
writing a simple product state is that quantum correla- 
tions can be built in easily by adding extra terms to Hi . 
We define 

IV^H ]>>>), (C4) 

where 

H = Hi + H 2 + H 3 + . . . , (C5) 

and i?2, H3, ■ ■ • denote many-body Jastrow factors. The 
correlated wave function (|C4|) is easy to use in VMC, as 
long as H is diagonal in the Ising basis \a). In this paper, 
we only consider two-body correlation terms, 

H 2 = - J2W S » S *j) + l{SziS zj ) 2 } . (C6) 

In principle, in Eq. (|C6j) . the sum can go over further- 
neighbor lattice sites, and the variational parameters (3 
and 7 may depend on the distance between sites. How- 
ever, inclusion of further-neighbor correlations are ex- 
pected to have a small effect on the ground state energy^ 
Because of this, and also, in order to have a number of 
variational parameters that is similar to the number of 
parameters used for the spin liquid wave functions, we 
consider only nearest-neighbor Jastrow factors here. 

The VMC algorithm can now be applied to Huse-Elser 
wave functions as outlined in Appendix The wave 
function is given by 

(c#) = e^ (Q) , (C7) 

where H(a) = (a\H\a). The Metropolis acceptance 
probability (|A2[) and the expectation values (|A3[) are 
straightforward to calculate. In contrast to the case of 
fcrmionic wave functions, no determinants or Pfaffians 
need to be evaluated or updated for this. 

In contrast to similar wave functions for spin S = 1/2, 
an important subtlety arises here in the generation of the 
random walk. For S = 1/2 and Sl ot = 0, the configu- 
rations a are restricted to states with an equal number 
of up and down spins. Therefore, the only admissible 
local Monte Carlo move is an exchange of two opposite 
spins. For S = 1, due to presence of the nematic state |0) 
with S z = 0, more local moves are possible. The Hilbert 
space for 5 = 1 and Sl ot = can be written as a direct 
sum of orthogonal subspaces (iVi-sectors) with a fixed 
number Ni = . . . N/2 of sites in configuration |1). The 
dimension of each iVi-sector is D(Ni) = ( N ) ( N ^ 1 ) ■ 
There exist two types of local moves in a random walk 
through the Ising configurations: those leaving Ni intact 
and those changing N± and moving to a different Ni- 
sector. The algorithm generating the random walk has 
to be unbiased with respect to moves between different 
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sectors such that each /Vi-sector is visited with proba- 
bility p(JV x ) = DiN^/J^nllDin). We have checked 
that such a distribution is accurately generated by the 
following procedure. We pick two sites at random and, 
depending on the states found on the sites, perform the 
following move: 

(i) |0)|1) or |0)|T): exchange the states. 

(ii) |1)|1): exchange states or change to |0)|0), each with 
probabilities 1/2. 

(hi) |0)|0): change state to |1)|1). 

(iv) |1)|1) or |1)|T): pick two different sites that are oc- 
cupied by unequal flavors and exchange them. 

In (iv), when the configurations |1)|1) or |1)|1) are en- 
countered, it is important to find two flavors to exchange, 
thus not changing iVi-sector. For example, if our algo- 
rithm rejected this case, and retried with another pair of 
sites, the random walk would be biased with respect to 
the distribution p(Ni), resulting in a higher probability 
for visiting sectors with smaller N\. 

Appendix D: Symmetries of KD- and 5(7(3)-models 

In this appendix, we elaborate on the symmetry prop- 
erties of the bilinear-biquadratic model (|14p and of the 
SU(3) ring-exchange model (l2Tj) investigated in this pa- 
per. 

Let us first discuss the ST/ (3) symmetry of these mod- 
els. Writing the Heisenberg exchange operator for spin 
S = 1, Eq. ([20]). in terms of the operators / = {f x ,f y ,f z ), 
we have 

V,j = Si • S, + (Si ■ S,r - 1 

— fiifbifbjfaj — fi ' (/» ' fj)fj ■ 

ab 

In this notation it is clear that Vij is invariant under a 
transformation / i— > Af where A is a general 3x3 unitary 
matrix. However, as discussed previously, the transfor- 
mation f a i-» e^f a with the same phase for all flavors 
does not change the corresponding spin state. Therefore, 
the relevant spin symmetry is SU(3) — U(3)/U(l), and 
we can take A € St/ (3). Similar to the operators f a that 
create these states, the spin states |o) transform in the 
fundamental representation of the SU (3) symmetry, by 
matrix multiplication with A. To find the action of the 
symmetry on spin operators, let us define 

ab 

where A p = (A" t b ), fj, = 1 ... 8, are the Gell-Mann matri- 
ces, generators of SU(3). Using 

[<W»]=I>?A> (D3) 

b 



it is clear that 

A f = e iad(Q) f = e iQ f e -iQ t (D4) 

for A = cxp{i J2,j, a^X^} and Q = a^Q^- Therefore, 
the spin operators 

S = -if A / (D5) 

transform as 

S i — ^ e^Se~^ (D6) 

under an SU (3) symmetry transformation. 

Rather than explicitly writing down all eight genera- 
tors of the SU(3) symmetry, Eq. (|D2|) . in spin language 
using the Gell-Mann basis, let us mention an equivalent 
set of generators. This set consist of the three spin rota- 
tion generators S a and the five independent quadrupolar 
operators Q ab = (S a S b + S b S a - 3)/2& 

The ring-exchange model, Eq. (|2ip. can be written in 
terms of Heisenberg exchange operators Vij. Therefore, 
it has the large SU(3) symmetry discussed above for all 
values of parameter a. The ifD-model, Eq. ([H]). enjoys 
the SU(3) symmetry only at the special point K = 1 
and D = in parameter space (where it is equivalent to 
the ring-exchange model at a — 0). Moving away from 
this special point, for general K but keeping D = 0, the 
symmetry is reduced to SO(3) spin rotation symmetry, 
generated by S. Finally, for D ^ 0, this symmetry is 
further reduced to U(l) spin rotation about the z-axis. 

When we move away from the SU (3) symmetric point 
along the line K = 1 and D / 0, the symmetry is re- 
duced to SU(2) on that line. Clearly, the spin rotation 
symmetry is reduced to S z as D ^ 0. To find the remain- 
ing unbroken generators, we need to determine the SU (3) 
generators that commute with the biquadratic term S|. 
These generators are S x S y + S y S x and S% — Sy. Hence, 
{S z , S x S y + S y S x , Si — Sy } are the three generators of 
an SU (2) symmetry of the model (fT4|) on the line K = 1. 

Let us briefly discuss the symmetry reasons behind 
the degeneracy of the correlated AFM and the nematic 
states, (fT5")) and ([T5|. on the line K = 1. In terms 
of spinon operators, the relevant symmetry generator is 
written as 

Si -Sl = flf x - flfy . (D7) 

From (ID4[) , we see that x and y states simply acquire an 
opposite phase under this transformation: f x t— > e llp f x , 
f y i y e~ %v fy. It is easy to check that the magnetic state 
(fT5"j) is mapped to the spin-nematic state (|T5|) for ip = 
7r/2, i.e., when f x i-> if x and f v i— > —if y . Furthermore, 
it is clear that the hopping term in ([9J and the Jastrow 
factors in (|10p are invariant under this transformation. 
Hence, the correlated ordered states are exactly mapped 
into each other by this transformation. 
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